//
// Solver.cpp for  in /home/herpec_j/Dropbox/Projets/raytracer-2016/Sources/Maths/
//
// Made by Jonathan
// Login   <ledey_l@epitech.net>
//
// Started on  Sun Sep 22 00:02:03 2013 Jonathan
// Last update Sun Sep 22 16:26:27 2013 Jonathan
//

#include	"Maths/Solver.hpp"

namespace	Maths
{
  int		Solver::SolveQuadric(long double rCoeff[3], long double rSolutions[2])
  {
    long double	p = rCoeff[1] / (2 * rCoeff[2]);
    long double	q = rCoeff[0] / rCoeff[2];
    long double	d = p * p - q;

    if (IS_ZERO(d))
      {
	rSolutions[0] = -p;
	return (1);
      }
    else if (d > 0)
      {
	d = sqrt(d);
	rSolutions[0] = d - p;
	rSolutions[1] = -d - p;
	return (2);
      }
    else
      {
	return (0);
      }
  }

  int		Solver::SolveCubic(long double rCoeff[4], long double rSolutions[3])
  {
    long double	a = rCoeff[2] / rCoeff[3];
    long double	b = rCoeff[1] / rCoeff[3];
    long double	c = rCoeff[0] / rCoeff[3];
    long double	a2 = a * a;
    long double	p = 1.0f / 3.0f * (-1.0f / 3.0f * a2 + b);
    long double	q = 1.0f / 2.0f * (2.0f / 27.0f * a * a2 - 1.0f / 3.0f * a * b + c);
    long double	p3 = p * p * p;
    long double	d = q * q + p3;
    int		solution;

    if (IS_ZERO(d))
      {
	if (IS_ZERO(q))
	  {
	    rSolutions[0] = 0;
	    solution = 1;
	  }
	else
	  {
	    long double	u = CBRT(-q);

	    rSolutions[0] = 2 * u;
	    rSolutions[1] = -u;
	    solution = 2;
	  }
      }
    else if (d < 0)
      {
	long double	phi = 1.0f / 3.0f * acos(-q / sqrt(-p3));
	long double	t = 2 * sqrt(-p);

	rSolutions[0] = t * cos(phi);
	rSolutions[1] = -t * cos(phi + PI / 3.0f);
	rSolutions[2] = -t * cos(phi - PI / 3.0f);
	solution = 3;
      }
    else
      {
	long double	sqrt_d = sqrt(d);
	long double	u = CBRT(sqrt_d - q);
	long double	v = -CBRT(sqrt_d + q);

	rSolutions[0] = u + v;
	solution = 1;
      }
    long double	sub = 1.0f / 3.0f * a;

    for (unsigned short index = 0;
	 index < solution;
	 ++index)
      {
	rSolutions[index] -= sub;
      }
    return (solution);
  }

  int		SolveQuartic(long double rCoeff[5], long double rSolutions[4])
  {
    long double	a = rCoeff[3] / rCoeff[4];
    long double	b = rCoeff[2] / rCoeff[4];
    long double	c = rCoeff[1] / rCoeff[4];
    long double	d = rCoeff[0] / rCoeff[4];
    long double	a2 = a * a;
    long double	p = -3.0f / 8.0f * a2 + b;
    long double	q = 1.0f / 8.0f * a2 * a - 1.0f / 2.0f * a * b + c;
    long double	r = -3.0f / 256.0f * a2 * a2 + 1.0f / 16.0f * a2 * b - 1.0f / 4.0f * a * c + d;
    long double	coeffs[4];
    int		solution;

    if (IS_ZERO(r))
      {
	coeffs[0] = q;
	coeffs[1] = p;
	coeffs[2] = 0;
	coeffs[3] = 1;
	solution = Solver::SolveCubic(coeffs, rSolutions);
	rSolutions[solution] = 0;
	++solution;
      }
    else
      {
	coeffs[0] = 1.0f / 2.0f * r * p - 1.0f / 8.0f * q * q;
	coeffs[1] = -r;
	coeffs[2] = -1.0f / 2.0f * p;
	coeffs[3] = 1;
	Solver::SolveCubic(coeffs, rSolutions);

	long double	z = rSolutions[0];
	long double	u = z * z - r;
	long double	v = 2 * z - p;

	if (IS_ZERO(u))
	  {
	    u = 0;
	  }
	else if (u > 0)
	  {
	    u = sqrt(u);
	  }
	else
	  {
	    return (0);
	  }
	if (IS_ZERO(v))
	  {
	    v = 0;
	  }
	else if (v > 0)
	  {
	    v = sqrt(v);
	  }
	else
	  {
	    return (0);
	  }
	coeffs[0] = z - u;
	coeffs[1] = q < 0 ? -v : v;
	coeffs[2] = 1;
	solution = Solver::SolveQuadric(coeffs, rSolutions);
	coeffs[0] = z + u;
	coeffs[1] = q < 0 ? v : -v;
	coeffs[2] = 1;
	solution += Solver::SolveQuadric(coeffs, rSolutions + solution);
      }
    long double	sub = 1.0f / 4.0f * a;

    for (unsigned short index = 0;
	 index < solution;
	 ++index)
      {
	rSolutions[index] -= sub;
      }
    return (solution);
  }
};
